Pseudotime analysis of T-cells in developing thymus

library(Seurat)
library(conos)
library(ggpubr)
library(tidyverse)
library(SingleCellExperiment)
library(chromVAR)
library(motifmatchr)
library(BSgenome.Hsapiens.UCSC.hg38)
# library(monocle3)
source("~/multiOmic_benchmark/utils.R")
source("~/multiOmic_benchmark/integrateBenchmark.R")
source("~/multiOmic_benchmark/preprocess/selectFeatures.R")

Based on the results of my benchmark, I set out to align expression and accessibility profiles from the F74 developing thymus dataset to detect changes in accessibility along pseudotime trajectories. While the benchmark was based on the task of label propagation, I here use the two most faithful methods (Seurat CCA and Conos) to achieve a common embedding of ATAC-seq and RNA-seq cells.

Load datasets.

seu.rna@assays$RNA@data

33694 x 8321 sparse Matrix of class "dgCMatrix"

   [[ suppressing 68 column names ‘AAACCTGAGTTCGATC_1’, ‘AAACCTGCAAGTTGTC_1’, ‘AAACCTGCAATACGCT_1’ ... ]]
   [[ suppressing 68 column names ‘AAACCTGAGTTCGATC_1’, ‘AAACCTGCAAGTTGTC_1’, ‘AAACCTGCAATACGCT_1’ ... ]]

                                                                                                                                                            
RP11-34P13.3  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ......
FAM138A       . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ......
OR4F5         . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ......
RP11-34P13.7  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ......
RP11-34P13.8  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ......
RP11-34P13.14 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ......
RP11-34P13.9  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ......

 ..............................
 ........suppressing 8253 columns and 33680 rows in show(); maybe adjust 'options(max.print= *, width = *)'
 ..............................

   [[ suppressing 68 column names ‘AAACCTGAGTTCGATC_1’, ‘AAACCTGCAAGTTGTC_1’, ‘AAACCTGCAATACGCT_1’ ... ]]

                                                                                                                                                         
AC023491.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ......
AC004556.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ......
AC233755.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ......
AC233755.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ......
AC240274.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ......
AC213203.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ......
FAM231B    . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ......

Filter genes with zero variance

rna.sce <- rna.sce[which(rna.gene.var > 0),]
atac.sce <- atac.sce[which(atac.gene.var > 0),]

rna.sce; atac.sce

class: SingleCellExperiment 
dim: 24510 8321 
metadata(0):
assays(3): counts cpm logcounts
rownames(24510): RP11-34P13.3 RP11-34P13.7 ... AC233755.1 AC240274.1
rowData names(0):
colnames(8321): AAACCTGAGTTCGATC_1 AAACCTGCAAGTTGTC_1 ... TTTGTCAAGCTGAACG_2 TTTGTCAGTATTAGCC_2
colData names(1): annotation
reducedDimNames(0):
spikeNames(0):
class: SingleCellExperiment 
dim: 31122 5793 
metadata(0):
assays(3): counts cpm logcounts
rownames(31122): A1BG A1BG-AS1 ... ZYX ZZEF1
rowData names(0):
colnames(5793): AAACGAAAGTGAACCG-1 AAACGAACATCGGCCA-1 ... TTTGTGTTCGATCGCG-1 TTTGTGTTCTGAGTAC-1
colData names(28): orig.ident nCount_ATAC ... nFeature_ACTIVITY ident
reducedDimNames(2): LSI UMAP
spikeNames(0):

Integration of T cells clusters

I re-run the integration based on the T cell subset. To select cells from the scATAC dataset, I take the SnapATAC clusters that best correspond to T-cells, based on label transfer.

tcells.sce.atac <- atac.sce[,which(as.numeric(atac.sce$seurat_clusters) %in% c(1:9))]

tcells.rna.ix <- which(rna.sce$annotation %in% c("DN","DP (Q)", "DP (P)", "SP (1)", "SP (2)"))
tcells.sce.rna <- rna.sce[,tcells.rna.ix]

tcells.sce.list <- list(RNA=tcells.sce.rna, ATAC=tcells.sce.atac)

## Make color palette 4 cell types
cell.types <- as.character(unique(tcells.sce.rna$annotation))
cell.type.pal <- brewer.pal(length(cell.types), "Set1") %>% rev() %>% setNames(cell.types)

Next, I select genes on which to perform integration. I take the union of the most variable features in the RNA dataset and the most covered features in the ATAC dataset

Remove cell cycle genes, that might interfere with pseudotime ordering

cell_cycle_genes <- read.table("~/annotations/cell_cycle_genes.tsv")$V1

integrate_features_union <- union(hvg.rna, hcg.atac)
integrate_features_union <- setdiff(integrate_features_union, cell_cycle_genes) 

## Select features in both datasets
integrate_features_union <- intersect(integrate_features_union, intersect(rownames(tcells.sce.list$ATAC), rownames(tcells.sce.list$RNA))) 

Visualize T cells in RNA dataset

tcells.seu.list <- map(tcells.sce.list, ~ as.Seurat(.x))

All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to LSI_All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP_

tcells.RNA.union <- tcells.seu.list$RNA
VariableFeatures(tcells.RNA.union) <- integrate_features_union
tcells.RNA.union <- ScaleData(tcells.RNA.union) %>% RunPCA() %>% RunUMAP(dims=1:40)

Centering and scaling data matrix

  |                                                                                  
  |                                                                            |   0%
  |                                                                                  
  |===============                                                             |  20%
  |                                                                                  
  |==============================                                              |  40%
  |                                                                                  
  |==============================================                              |  60%
  |                                                                                  
  |=============================================================               |  80%
  |                                                                                  
  |============================================================================| 100%
The following 176 features requested have zero variance (running reduction without them): TM4SF18, MEOX2, ASGR2, DMRT2, KIR2DL4, CACNA1B, FREM2, MYZAP, TOX3, PRIMA1, PPP2R2B, PTPRZ1, CMTM5, TREM2, CALB2, TMTC1, CDH6, PKHD1L1, ISLR2, KRT7, MYH8, DCX, CCSER1, MSLN, CDH19, EYA4, P2RY12, PLA2G2D, ADAMTS16, KCNJ5, PAPPA, QRFPR, ADGRG2, NFIA-AS2, LUZP2, NALCN, FABP7, TSPEAR-AS1, GPR1, CD300LF, A4GALT, AP000439.1, CEACAM3, SLITRK2, KIR3DL1, AC147651.1, MUC15, FFAR4, C11orf53, KLK11, WIF1, CFHR1, TMEM233, FOLR3, GRID2, GALNT15, FAM19A1, TAC1, PTCHD4, CCL11, ZNF385B, GABRA1, ADH1B, LINC01048, LMOD1, TNNT1, ACTL6B, SEMA3E, HOXD-AS2, NTF4, SLC35F1, MANEAL, FEV, RAB3C, SYT9, HMCN2, ZFHX4-AS1, PCSK9, MMP12, ABCB11, AC002066.1, ADAM23, ADCY5, ADGB, ADGRB1, ANKRD29, ATP2B2, C4orf45, CCDC33, CDH18, CDK15, CERS1, CFAP161, CLRN1-AS1, CNBD1, CNGB3, CNTN5, COL8A2, CTNND2, DCC, DSCAM, EPHA6, FGF12, FSTL5, GABBR2, GABRB1, GABRG3, GALNT13, GHR, GLIS3, GNG12-AS1, GNGT1, GRIN2B, GRM1, GRM8, HMCN1, HPSE2, HS3ST5, HS6ST3, IL1RAPL2, ITGBL1, KCNB1, KCNH7, LINC00639, LINC01169, LINC01170, LINC01182, LINC01317, LINC01505, LINGO1, LINGO2, LRRC4C, LRRTM4, LY75-CD302, MAP6, MLIP, MYRIP, NEBL, NKAIN2, NKAIN3, NRG1, OPCML, PCDH15, PCDHGA4, PDZRN4, PIK3C2G, PKNOX2, PNPLA1, POU6F2, PPFIA2, PRR5-ARHGAP8, PTPN5, PTPRT, RALYL, RERG, RGS7, RIMS1, RIN2, RYR3, SHANK2, SNTG1, SORCS3, STON1-GTF2A1L, SYN2, TENM2, THSD4, TRPC6, TSPEAR, TVP23C-CDRT4, UNC80, USH1C, VAX2, VWA3B, XKR4, ZBTB7C, ZNF804BPC_ 1 
Positive:  SATB1, PTPRC, MTSS1, APBB1IP, LYST, SH2D1A, CAMK4, LBH, FBLN5, LEF1 
       PLEKHG1, VOPP1, CD247, MXD1, TCF12, ARHGEF7, ALDH1A2, NFATC3, TBC1D19, SYTL3 
       ADAMTS17, ANO6, DAPK1, CALN1, THEMIS, PITPNM2, NLGN4X, MAPRE2, GALNT7, ZNF280D 
Negative:  ENO1, GSTP1, FABP5, TMSB10, SMS, PKM, NME4, VIM, RPL37A, YWHAQ 
       NCL, IGFBP2, CAPG, NDUFA12, TRDC, PGK1, PARVB, LDHB, ATOX1, SELL 
       CDC123, NUDC, IGLL1, UBE2N, FXYD2, GMPS, C20orf27, SLC25A39, ANXA1, C12orf75 
PC_ 2 
Positive:  RPL37A, EIF3H, LDHB, TBCA, IL32, SERGEF, SMPD3, ITGAE, AATF, FBLN5 
       NDUFA12, SOD1, DNAJC15, C12orf75, MRPL33, TMSB10, HNRNPC, CCDC57, SMCO4, GDI2 
       OLA1, ALDH1A2, COX7A2L, CST3, CYSTM1, ATOX1, GYPC, PDCD6, FABP5, RALY 
Negative:  MBNL1, ITPR2, PTPRC, MTRNR2L12, ADAM10, MSI2, CDK6, SELL, JCHAIN, RNF213 
       MME, TCF12, BPTF, BCL11B, NIPBL, MACF1, PIK3R1, HIVEP3, SOCS2, GALNT7 
       IKZF2, RUNX1, SPTBN1, PRRC2B, BCL2, DIAPH1, MYCBP2, SLC38A1, MBP, RUFY3 
PC_ 3 
Positive:  HLA-B, TOX2, COTL1, CTSW, KLRB1, CD40LG, SIRPG, GZMM, CLDN1, CD74 
       ITM2A, GBP2, BACH2, HPGD, PDE4D, TNFRSF1B, S100A10, CLEC2D, XCL1, GFOD1 
       CRTAM, MATK, DENND2D, PDCD1, GIMAP4, STAT1, ZNF683, CD226, HLA-A, TNFRSF25 
Negative:  TFDP2, JCHAIN, ATP6AP1L, DEFA6, MME, GALNT2, PCGF5, ADGRG1, GLIPR1, MSI2 
       NINL, CEP70, CDK6, PTPN2, FXYD2, TRDC, GSTP1, SELL, SOCS2, RGPD3 
       SMPD3, FABP5, LYST, SSBP2, PITPNM2, DLEU7, UBE2E1, NUCB2, PPP1R1C, BCL11A 
PC_ 4 
Positive:  SMPD3, C12orf75, TUBA1C, LCP1, HIST1H2AB, SMS, SMCO4, GMPS, FABP5, XPO1 
       IGFBP2, CCND3, TAF15, RGS3, SREBF2, YWHAQ, TMSB15A, ABCD3, NCL, PHIP 
       GNAS, EPB41L2, SYNE2, STAG2, CNTLN, VIM, NUP210, MCTP1, NUDC, NFATC3 
Negative:  JCHAIN, SELL, DEFA6, SOCS2, GLIPR1, ATP6AP1L, MME, RGPD3, XG, DPP4 
       IFI6, ENAM, GNAL, EVL, NEGR1, FRMPD2, EVA1A, PIK3CD, GALNT2, NDST3 
       MBP, RNF144A, FXYD2, COL1A2, SMIM24, NDFIP2, OXNAD1, ACTN1, LSP1, LINC00861 
PC_ 5 
Positive:  HPGD, BACH2, TOX2, ITM2A, GZMM, ST6GAL1, ARAP2, CD96, TGFBR2, LZTFL1 
       VIM, EVL, CD44, SATB1, IL2, TUSC3, CYTIP, PDE4D, MAD1L1, SLC2A3 
       ARHGAP31, GPR183, PGK1, ANKRD44, STK4, BCL2, GNG2, IKZF1, ICOS, ETS1 
Negative:  XCL1, TNFRSF9, CTSW, XCL2, GBP2, TNFRSF1B, S100A4, GNG4, NPW, CD74 
       HOPX, IGFBP4, CD151, SH3BGRL2, NFATC1, LYST, ATP9A, NR4A2, LINC01480, CLEC2D 
       MIR3142HG, LINC01281, PREX1, HDAC9, PHACTR1, MAP7, CST3, CXCR3, PITPNM2, PDCD1 
09:42:53 UMAP embedding parameters a = 0.9922 b = 1.112
09:42:53 Read 7101 rows and found 40 numeric columns
09:42:53 Using Annoy for neighbor search, n_neighbors = 30
09:42:53 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
09:42:55 Writing NN index file to temp file /tmp/RtmpInlUFP/file4d1e59aaa998
09:42:55 Searching Annoy index using 1 thread, search_k = 3000
09:42:57 Annoy recall = 100%
09:42:59 Commencing smooth kNN distance calibration using 1 thread
09:43:03 Initializing from normalized Laplacian + noise
09:43:03 Commencing optimization for 500 epochs, with 314138 positive edges
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
09:43:22 Optimization finished

DimPlot(tcells.RNA.union, group.by = "annotation", label=TRUE) + ggtitle("RNA - feature union")

Visualize markers

t.cell.markers <- list(known.markers = c("CD34", "IGLL1", "TRGC2", "TRDC", "PTCRA", "TRBC2", "TRAC", "CD4", "CD8A", "CD8B"),
                       chemokine.receptors = c("CCR9", "CCR7"),
                       tcr.activation = c("CD5", "CD27"),
                       proliferation=c("PCNA", "CDK1", "MKI67"),
                       cyclin.D = c("CCND2", "CCND3"),
                       recombination=c("RAG1", "RAG2"),
                       apoptosis=c("HRK","BMF", "TP53INP1"),
                       stage.markers = c("ST18", "HIVEP3", "RGPD3", "SMPD3", "AQP3", "RORC", "SATB1", "TOX2")
                       ) 
# FeaturePlot(tcells.RNA.ref, features = t.cell.markers$known.markers, cols = viridis::viridis(n=10))
FeaturePlot(tcells.RNA.union, features = t.cell.markers$known.markers, cols = viridis::viridis(n=10))

Visualize T cells in ATAC dataset: is the trajectory visible in the binary matrix?

tcells.ATAC.union <- tcells.seu.list$ATAC
# tcells.ATAC.union <- NormalizeData(tcells.ATAC.union)
VariableFeatures(tcells.ATAC.union) <- integrate_features_union
tcells.ATAC.union <- RunLSI(tcells.ATAC.union, n=50, scale.max = NULL)

RunLSI is being moved to Signac. Equivalent functionality can be achieved via the Signac::RunTFIDF and Signac::RunSVD functions; for more information on Signac, please see https://github.com/timoast/SignacRunLSI is being moved to Signac. Equivalent functionality can be achieved via the Signac::RunTFIDF and Signac::RunSVD functions; for more information on Signac, please see https://github.com/timoast/SignacRunLSI is being moved to Signac. Equivalent functionality can be achieved via the Signac::RunTFIDF and Signac::RunSVD functions; for more information on Signac, please see https://github.com/timoast/SignacPerforming TF-IDF normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Running SVD on TF-IDF matrix
Scaling cell embeddings
Cannot add objects with duplicate keys (offending key: LSI_), setting key to 'lsi_'

tcells.ATAC.union <- RunUMAP(tcells.ATAC.union, reduction = "lsi", dims = 1:50)

09:43:59 UMAP embedding parameters a = 0.9922 b = 1.112
09:43:59 Read 4977 rows and found 50 numeric columns
09:43:59 Using Annoy for neighbor search, n_neighbors = 30
09:43:59 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
09:44:00 Writing NN index file to temp file /tmp/RtmpInlUFP/file4d1e64c0f0d8
09:44:00 Searching Annoy index using 1 thread, search_k = 3000
09:44:02 Annoy recall = 100%
09:44:04 Commencing smooth kNN distance calibration using 1 thread
09:44:08 Initializing from normalized Laplacian + noise
09:44:08 Commencing optimization for 500 epochs, with 172060 positive edges
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
09:44:21 Optimization finished
Cannot add objects with duplicate keys (offending key: UMAP_), setting key to 'umap_'

DimPlot(tcells.ATAC.union, reduction = "umap", group.by = "seurat_clusters", label = TRUE) + ggtitle("ATAC gmat")

Run CCA and Conos

Run CCA

sce.list <- tcells.sce.list
reference = "RNA"
query = "ATAC" 
seurat.list <- imap(sce.list, ~ as.Seurat(.x, assay=.y))

All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to LSI_All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP_

seurat.list <- imap(seurat.list, ~ RenameCells(.x, add.cell.id=.y))
## Scale data
seurat.list <- map(seurat.list, ~ ScaleData(.x))

Centering and scaling data matrix

  |                                                                                                                                                     
  |                                                                                                                                               |   0%
  |                                                                                                                                                     
  |======                                                                                                                                         |   4%
  |                                                                                                                                                     
  |===========                                                                                                                                    |   8%
  |                                                                                                                                                     
  |=================                                                                                                                              |  12%
  |                                                                                                                                                     
  |=======================                                                                                                                        |  16%
  |                                                                                                                                                     
  |=============================                                                                                                                  |  20%
  |                                                                                                                                                     
  |==================================                                                                                                             |  24%
  |                                                                                                                                                     
  |========================================                                                                                                       |  28%
  |                                                                                                                                                     
  |==============================================                                                                                                 |  32%
  |                                                                                                                                                     
  |===================================================                                                                                            |  36%
  |                                                                                                                                                     
  |=========================================================                                                                                      |  40%
  |                                                                                                                                                     
  |===============================================================                                                                                |  44%
  |                                                                                                                                                     
  |=====================================================================                                                                          |  48%
  |                                                                                                                                                     
  |==========================================================================                                                                     |  52%
  |                                                                                                                                                     
  |================================================================================                                                               |  56%
  |                                                                                                                                                     
  |======================================================================================                                                         |  60%
  |                                                                                                                                                     
  |============================================================================================                                                   |  64%
  |                                                                                                                                                     
  |=================================================================================================                                              |  68%
  |                                                                                                                                                     
  |=======================================================================================================                                        |  72%
  |                                                                                                                                                     
  |=============================================================================================================                                  |  76%
  |                                                                                                                                                     
  |==================================================================================================================                             |  80%
  |                                                                                                                                                     
  |========================================================================================================================                       |  84%
  |                                                                                                                                                     
  |==============================================================================================================================                 |  88%
  |                                                                                                                                                     
  |====================================================================================================================================           |  92%
  |                                                                                                                                                     
  |=========================================================================================================================================      |  96%
  |                                                                                                                                                     
  |===============================================================================================================================================| 100%
Centering and scaling data matrix

  |                                                                                                                                                     
  |                                                                                                                                               |   0%
  |                                                                                                                                                     
  |====                                                                                                                                           |   3%
  |                                                                                                                                                     
  |=========                                                                                                                                      |   6%
  |                                                                                                                                                     
  |=============                                                                                                                                  |   9%
  |                                                                                                                                                     
  |==================                                                                                                                             |  12%
  |                                                                                                                                                     
  |======================                                                                                                                         |  16%
  |                                                                                                                                                     
  |===========================                                                                                                                    |  19%
  |                                                                                                                                                     
  |===============================                                                                                                                |  22%
  |                                                                                                                                                     
  |====================================                                                                                                           |  25%
  |                                                                                                                                                     
  |========================================                                                                                                       |  28%
  |                                                                                                                                                     
  |=============================================                                                                                                  |  31%
  |                                                                                                                                                     
  |=================================================                                                                                              |  34%
  |                                                                                                                                                     
  |======================================================                                                                                         |  38%
  |                                                                                                                                                     
  |==========================================================                                                                                     |  41%
  |                                                                                                                                                     
  |===============================================================                                                                                |  44%
  |                                                                                                                                                     
  |===================================================================                                                                            |  47%
  |                                                                                                                                                     
  |========================================================================                                                                       |  50%
  |                                                                                                                                                     
  |============================================================================                                                                   |  53%
  |                                                                                                                                                     
  |================================================================================                                                               |  56%
  |                                                                                                                                                     
  |=====================================================================================                                                          |  59%
  |                                                                                                                                                     
  |=========================================================================================                                                      |  62%
  |                                                                                                                                                     
  |==============================================================================================                                                 |  66%
  |                                                                                                                                                     
  |==================================================================================================                                             |  69%
  |                                                                                                                                                     
  |=======================================================================================================                                        |  72%
  |                                                                                                                                                     
  |===========================================================================================================                                    |  75%
  |                                                                                                                                                     
  |================================================================================================================                               |  78%
  |                                                                                                                                                     
  |====================================================================================================================                           |  81%
  |                                                                                                                                                     
  |=========================================================================================================================                      |  84%
  |                                                                                                                                                     
  |=============================================================================================================================                  |  88%
  |                                                                                                                                                     
  |==================================================================================================================================             |  91%
  |                                                                                                                                                     
  |======================================================================================================================================         |  94%
  |                                                                                                                                                     
  |===========================================================================================================================================    |  97%
  |                                                                                                                                                     
  |===============================================================================================================================================| 100%

## Calculate CCA anchors
transfer.anchors <- FindTransferAnchors(reference = seurat.list[[reference]], 
                                        query = seurat.list[[query]],
                                        features = integrate_features_union, 
                                        reduction = "cca")

Running CCA on different assaysRunning CCA
Merging objects
Finding neighborhoods
Finding anchors
    Found 19433 anchors
Filtering anchors
    Retained 1616 anchors
Extracting within-dataset neighbors

## Impute expression profiles for ATAC cells (for all genes, not just integration features)
refdata <- GetAssayData(seurat.list$RNA, assay = "RNA", slot = "data")
imputation <- TransferData(anchorset = transfer.anchors, refdata = refdata, weight.reduction = seurat.list$ATAC[["LSI"]])

Finding integration vectors
Finding integration vector weights
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Transfering 24510 features onto reference data

## Merge datasets and co-embed
seurat.list$ATAC[["RNA"]] <- imputation
coembed <- merge(x = seurat.list$RNA, y = seurat.list$ATAC)

coembed <- ScaleData(coembed, features = integrate_features_union, do.scale = FALSE)

Centering data matrix

  |                                                                                                                                                     
  |                                                                                                                                               |   0%
  |                                                                                                                                                     
  |=============================                                                                                                                  |  20%
  |                                                                                                                                                     
  |=========================================================                                                                                      |  40%
  |                                                                                                                                                     
  |======================================================================================                                                         |  60%
  |                                                                                                                                                     
  |==================================================================================================================                             |  80%
  |                                                                                                                                                     
  |===============================================================================================================================================| 100%

coembed <- RunPCA(coembed, features = integrate_features_union, verbose = FALSE)
coembed <- RunUMAP(coembed, dims = 1:30)

09:52:05 UMAP embedding parameters a = 0.9922 b = 1.112
09:52:05 Read 12078 rows and found 30 numeric columns
09:52:05 Using Annoy for neighbor search, n_neighbors = 30
09:52:05 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
09:52:07 Writing NN index file to temp file /tmp/RtmpInlUFP/file4d1e3a771214
09:52:07 Searching Annoy index using 1 thread, search_k = 3000
09:52:12 Annoy recall = 100%
09:52:15 Commencing smooth kNN distance calibration using 1 thread
09:52:19 Initializing from normalized Laplacian + noise
09:52:19 Commencing optimization for 200 epochs, with 559810 positive edges
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
09:52:33 Optimization finished

coembed <- AddMetaData(coembed, metadata = ifelse(colnames(coembed) %in% colnames(seurat.list[[reference]]), reference, query), col.name = "tech")

DimPlot(coembed, group.by = c('tech', "annotation"))

FeaturePlot(coembed, features = t.cell.markers$stage.markers, split.by = "tech", cols = viridis::viridis(n=100)) + ggtitle("Stage Markers")

FeaturePlot(coembed, features = t.cell.markers$known.markers, split.by = "tech", slot = "data", cols = viridis::viridis(n=100))

FeaturePlot(coembed, features = t.cell.markers$recombination, split.by = "tech", slot = "data", cols = viridis::viridis(n=100)) 

Transfer labels on ATAC dataset

celltype.predictions <- TransferData(anchorset = transfer.anchors, 
                                     refdata = seurat.list[[reference]]$annotation, 
                                     weight.reduction = seurat.list$ATAC[["LSI"]])

Finding integration vectors
Finding integration vector weights
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Predicting cell labels

coembed <- AddMetaData(coembed, metadata = celltype.predictions)
coembed@meta.data %<>%
  rownames_to_column() %>%
  dplyr::mutate(annotation=ifelse(is.na(predicted.id) , annotation, NA)) %>%
  column_to_rownames()

coembed@meta.data <-
  coembed@meta.data %>%
  rownames_to_column() %>%
  dplyr::mutate(annotation=ifelse(is.na(annotation) & prediction.score.max > 0.5, predicted.id, annotation)) %>%
  column_to_rownames()
CombinePlots(
  list(DimPlot(coembed, group.by = c("predicted.id"), cols = cell.type.pal) + ggtitle("prediction"),
  DimPlot(coembed, group.by = c("annotation"), cols = cell.type.pal) + ggtitle("Original + prediction")),
  legend = "top"
  )

FeaturePlot(coembed, features = "prediction.score.max", cells = which(coembed$tech=="ATAC")) + scale_color_viridis_c()

Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.

Run Pseudotime analysis

Identify cell of origin based on IGLL1 and CD34

FeaturePlot(coembed, features = c("IGLL1", "CD34"), split.by = "tech", slot = "data", cols = viridis::viridis(n=100))

cell.oo <-
  coembed@meta.data %>% 
  rownames_to_column("cell") %>%
  mutate(IGLL1=coembed@assays$RNA@counts["IGLL1",cell]) %>%
  select(cell, annotation, IGLL1) %>%
  arrange(-IGLL1) %>%
  filter(annotation=="DN") %>%
  top_n(1, IGLL1) %>%
  pull(cell)

coembed@reductions$umap@cell.embeddings %>%
  as.tibble(rownames="cell") %>%
  mutate(cell.oo = ifelse(cell %in% cell.oo, T, F)) %>%
  ggplot(aes(UMAP_1, UMAP_2)) +
  geom_point(color="grey50") +
  geom_point(data=. %>% filter(cell.oo),color='red') +
  ggrepel::geom_text_repel(data=. %>% filter(cell.oo), aes(label="cell of origin"), color='red') +
  theme_cowplot() 

`as.tibble()` is deprecated, use `as_tibble()` (but mind the new semantics).
This warning is displayed once per session.

coembed <- AddMetaData(coembed, ifelse(colnames(coembed)==cell.oo, TRUE, FALSE), col.name = "iroot_cell")

  

I infer pseudotime using the diffusion pseudotime algorithm as implemented in scanpy. Making an R/reticulate wrapper for this function would be nice, but for now, see multiOmic_benchmark/DPT_tcells.ipynb.

Read scanpy output and save in R object.

dpt <- read.csv('~/my_data/Tcells_CCA_integration_20191127_scanpy_dpt.csv') %>%
  select(X, dpt_pseudotime)

coembed <- AddMetaData(coembed, column_to_rownames(dpt, 'X'))
saveRDS(coembed, "~/my_data/Tcells_CCA_integration_seurat_20191127.Rmd")

coembed <- readRDS("~/my_data/Tcells_CCA_integration_seurat_20191127.Rmd")

Visualize pseudotime

Save figure

Check expression of markers along pseudotime

coembed@assays$RNA@data[t.cell.markers$known.markers, ] %>%
  as.matrix() %>%
  reshape2::melt(varnames=c("gene", "cell")) %>%
  group_by(gene) %>%
  # mutate(value=(value-mean(value))/sd(value)) %>%
  left_join(coembed@meta.data[,"dpt_pseudotime", drop=F] %>% rownames_to_column("cell")) %>%
  mutate(pseudotime.rank=dense_rank(dpt_pseudotime)) %>%
  ggplot(aes(pseudotime.rank, gene)) +
  geom_tile(aes(fill=value)) +
  scale_fill_viridis_c()

Bin pseudotime and visualize cell type composition

Correlation between global accessibility and pseudotime ordering.

Motif analysis

I initially wanted to call peaks from SnapATAC clusters, then build a cell x peak matrix on those detected peaks, but SnapATAC/MACS2 were drivin me nuts.

Alternative: load peak matrix from cellranger and add to snap object

snap.pmat <- createSnapFromPmat(mat=t(filt.peaks[,snap.out@barcode]), barcodes=snap.out@barcode, peaks=peaks.gr)
snap.pmat

number of barcodes: 5793
number of bins: 0
number of genes: 0
number of peaks: 93607
number of motifs: 0

Calculating deviations in TF accessibility using ChromVAR. This is a measure of how much is motif accessibility in each cell is enriched compared to all the cells and general cell coverage. While SnapATAC has an wrapper around ChromVAR that outputs the deviation matrix, I just take the code from that function and run every step separately to keep the useful outputs and statistics of chromVAR.

snap.pmat = makeBinary(snap.pmat, "pmat")

obj = snap.pmat
input.mat="pmat"
min.count=10
species="Homo sapiens"
genome=BSgenome.Hsapiens.UCSC.hg38

data.use = obj@pmat
peak.use = obj@peak

ncell = nrow(data.use)

idy = which(Matrix::colSums(data.use) >= min.count)
data.use = data.use[,idy,dropping=TRUE]
    
peak.use = peak.use[idy]

rse <- SummarizedExperiment(
        assays = list(counts = t(data.use)), 
                 rowRanges = peak.use, 
                 colData = DataFrame(Cell_Type=1:nrow(data.use), depth=Matrix::rowSums(data.use))
    );
rse <- addGCBias(rse, genome = genome);
motifs <- getJasparMotifs(collection = "CORE", species=species);
motif_mm <- matchMotifs(motifs, rse, genome = genome);
dev <- computeDeviations(object = rse, annotations = motif_mm);

'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading

var <- computeVariability(dev)

'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading

plotVariability(var, use_plotly = F, n = 8)

Save

rowData(dev) %<>%
  as.tibble(rownames="motif") %>%
  full_join(var) %>%
  column_to_rownames('motif') %>%
  DataFrame()

Joining, by = "name"
Column `name` joining character vector and factor, coercing into character vector

Visualize motif enrichment

Compare motif accessibility and TF gene expression

genex.df <- coembed@assays$RNA@data[tf.topvar[which(tf.topvar %in% rownames(coembed@assays$RNA))],] %>%
  as.matrix() %>%
  reshape2::melt(varnames=c("gene", 'cell')) %>%
  filter(str_detect(cell, "RNA_")) %>%
  left_join(dpt.df) %>%
  arrange(dpt_bin) %>%
  mutate(tech="RNA")

Joining, by = "cell"
Column `cell` joining factor and character vector, coercing into character vector

bind_rows(access.df, genex.df) %>%
  filter(gene %in% c("SPI1", "RUNX2","RUNX3", 'TCF7L2', "E2F4")) %>%
    filter(!str_detect(cell, "^ATAC_") | tech=="ATAC") %>%
  # group_by(tech, gene) %>%
  # mutate(value=scale(value)) %>%
  # ungroup() %>%
  # drop_na() %>%
# filter(tech=="RNA") %>%
  drop_na(dpt_bin) %>%
  ggplot(aes(dpt_bin, value)) +
  geom_point(aes(color=annotation), alpha=0.2) +
  facet_grid(tech~gene, scales = "free_y") +
  geom_smooth()

binding character and factor vector, coercing into character vector

bind_rows(access.df, genex.df) %>%
  filter(gene %in% c("ELK3", "JUNB", "FOS")) %>%
  filter(!str_detect(cell, "^ATAC_") | tech=="ATAC") %>%
  # group_by(tech, gene) %>%
  # mutate(value=scale(value)) %>%
  # ungroup() %>%
  drop_na(dpt_bin) %>%
  ggplot(aes(dpt_bin, value)) +
  geom_point(aes(color=annotation), alpha=0.2) +
  facet_grid(tech~gene, scales = "free_y") +
  geom_smooth() +
  scale_color_manual(values=cell.type.pal)

binding character and factor vector, coercing into character vector

Pseudotime lag between DP(Q) in accessibility and gene expression

The DP (Q) cluster in the ATAC cells is scored with high confidence

DimPlot(cca.obj, group.by=c("annotation","tech"), reduction = "cca", dims = 1:2)
DimPlot(cca.obj, group.by=c("annotation","tech"), reduction = "cca", dims = 3:4)

DimPlot(cca.obj, group.by=c("annotation","tech"), reduction = "cca", dims = 5:6)

plot_grid(
  DimPlot(dpq.coembed, group.by = c("tech"), reduction = "pca") ,
  FeaturePlot(dpq.coembed, feature="dpt_pseudotime", reduction = "pca") + scale_color_viridis_c() 
  )

Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.

FeaturePlot(dpq.coembed, feature=unique(topPC1.df$label)[1:6], reduction = "umap") 

The following requested variables were not found: NA

FeaturePlot(dpq.coembed, feature="FABP5", reduction = "pca", slot = "raw.data") 

Error in slot(object = object, name = slot) : 
  no slot of name "raw.data" for this object of class "Assay"

---